Nonequilibrium transport through molecular junctions in the quantum regime 
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We consider a quantum dot, affected by a local vibrational mode and contacted to macroscopic 
leads, in the nonequilibrium steady-state regime. We apply a variational Lang-Firsov transforma- 
tion and solve the equations of motion of the Green functions in the Kadanoff-Baym formalism up 
to second order in the interaction coefficients. The variational determination of the transformation 
parameter through minimization of the thermodynamic potential allows us to calculate the elec- 
tron/polaron spectral function and conductance for adiabatic to anti-adiabatic phonon frequencies 
and weak to strong electron-phonon couplings. We investigate the qualitative impact of the quasi- 
particle renormalization on the inelastic electron tunneling spectroscopy signatures and discuss the 
possibility of a polaron induced negative differential conductance. In the high-voltage regime we 
find that the polaron level follows the lead chemical potential to enhance resonant transport. 

PACS numbers: 72.10.-d, 71.38.-k, 73. 21. La, 73.63.Kv 
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I. INTRODUCTION 

Recent advances in nano-technology have made pos- 
sible the creation of electronic devices with the active 
element being a single organic molecule. Such molecular 
junctions may be an alternative to semi-conductor tech- 
nology in the search for further miniaturization and novel 
transport properties. They can be described as quantum 
dots, i.e., as systems of finite size coupled to macroscopic 
leads acting as charge reservoirs. As with metallic or 
semi-conducting junctions, energy level quantization de- 
termines transport. In addition, when being occupied by 
charge carriers, molecular quantum dots are susceptible 
to structural changes that may be induced by the interac- 
tion with optical phonons. As a consequence, vibrational 
signatures show up in the current- voltage characteristics. 
Moreover, they render inelastic tunneling spectroscopy 
(lETS) the primary experimental tool for the identifica- 
tion and characterization of molecular quantum dots.^'^ 

For a thorough understanding of the underlying trans- 
port mechanisms suitable theoretical models have to be 
studied. The simplest one is based on a modified Fano- 
Anderson model where the static impurity is replaced by 
a single site coupled to a local phonon mode. Then the 
current is given by the interacting dot spectral function 
and the voltage bias between the noninteracting macro- 
scopic leads^? The transport properties of the system 
strongly depend on the relative time scales of the elec- 
tronic and phononic subsystems)^ 

In the regime of fast electron motion and weak 
electron-phonon (EP) coupling standard perturbation 
theory appliesj^^— Here, lETS signatures result from 
the interference of (quasi)elastic and inelastic tunneling 
processes)^ The calculated line shapes in the total cur- 
rent are found to be especially sensitive to changes in the 
dot-lead coupling parameter and the dot level energy.— 
In general, both these quantities should be affected by 
conformational changes of the molecule. In the equilib- 



rium situation, the question remains whether vibrational 
coupling leads to a broadening^ or narrowing^ of the lin- 
ear conductance resonance as a function of the dot level. 

On the other hand, in molecular quantum dots the 
vibrational frequency can be larger than the kinetic en- 
ergy of incident electrons. From the study of the Hol- 
stein molecular crystal modelji^ it is well-known that 
in this regime, strong EP interaction may heavily re- 
duce the "mobility" of the electrons through the for- 
mation of small polarons (electrons dressed by phonon 
clouds) ji^ii^ Consequently, for quantum dots, the forma- 
tion of a local polaron is considered a possible mechanism 
for the observed nonlinear transport properties, such as 
hysteresis, negative differential conductance (NDC) and 
switching. ^^"^^ Approaches based on the application of 
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to the Hamiltonian 



suggest that the vibrational structure of the polaron 
state is revealed by distinct steps in the current-voltage 
signal Here, electron transport takes place via res- 
onant tunneling through phonon sidebands. 

In this paper, we investigate steady-state transport 
through molecular quantum dots for small-to-large dot- 
lead coupling and weak-to-strong EP interaction. Using 
the Meir-Wingreen current formula;^ our main task is the 
determination of the interacting electronic spectral func- 
tion of the quantum dot. As the background of our cal- 
culations we choose the formalism of Kadanoff-Baymj^ 
which relies on the correspondence of the nonequilibrium 
Green functions of complex times to the real-time re- 
sponse functions. Starting from the Dyson equation, the 
general steady-state equations for the response functions 
will be deduced. The solution of the latter equations 
will lead to a nonequilibrium spectral function which has 
a form analogous to the equilibrium one. The dot self- 
energy determining the spectral function will be calcu- 
lated from the equations of motion of the Green functions 
up to a second order in the interaction coefficients. 

Our approach is based on a variational Lang-Firsov 
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transformation, which was developed for Holstein po- 
larons at finite densities^ and recently applied to the 
molecular quantum dot in equilibriumj^S. We extend these 
calculations to the nonequilibrium situation and to fi- 
nite temperatures, whereby the dot self-energy will be 
calculated self-consistently to account for the density- 
dependent oscillator shift. The variational parameter 
of the Lang-Firsov transformation is determined numer- 
ically via the minimization of the thermodynamic po- 
tential. In this way we are able to interpolate between 
the self-consistent Born approximation (SCBA)iii2£ and 
the small-polaron approach^^ previously used in the weak 
and strong EP coupling limits. We note that already in 
the equilibrium case, our variational calculation intro- 
duces important corrections to the corresponding spec- 
tral functions, that determine the conductance in the lin- 
ear response theory. We reexamine the low-temperature 
equilibrium quantum dot system and analyze the occur- 
rence of high-temperature phonon sidebands in the linear 
conductance. 



In the nonequilibrium situation we show the impact of 
the optimal polaron state on the lETS signatures men- 
tioned above. For comparable electronic and phononic 
time scales, we study the crossover from coherent tunnel- 
ing to sequential hopping via a transient polaron state, 
where the interplay of both resonant and off-resonant 
multiphonon processes leads to complicated electron tun- 
neling spectra. Recently La Magna and Deretzis^^ ap- 
plied a similar variational ansatz to an effective electron 
Hamiltonian and found polaron-formation-induced NDC. 
Considering the dependence of the current-voltage char- 
acteristics on the full spectral function, we critically dis- 
cuss this effect. 



The paper is organized as follows: Sec. Ill Al intro- 
duces the model Hamiltonian and describes the varia- 
tional Lang-Firsov transformation. In Sees. Ill Bj and lll Cl 
a formal steady-state solution to the equations of motion 
is presented. In Sec. IIIDI we derive an approximation 
to the polaronic self-energy that is self-consistent and 
depends on the variational parameter. The latter is de- 
termined from the numerical minimization of the ther- 
modynamic potential that is deduced in Sec. Ill El Sec- 
tion IIIFI gives the relation between the electronic and 
polaronic spectral functions. In Sec. IIIGI the general 
current formula for arbitrary voltage is discussed and 
the special case of linear conductance is mentioned. Sec- 
tion |llT] presents our numerical results and Sec. |IV] sum- 
marizes. 



II. THEORY 
A. General equations 

Our considerations are based on the standard Hamil- 
tonian of the single-site quantum dot model, 

H = {A- ^i)dU - gujQdU{b^ + b) + Loob^b (1) 

k.a k.a 

Here, the quantum dot is represented by the energy level 
A, with the fermionic creation (destruction) operator d) 
[d). The dot is coupled to a local phonon mode 6^^' of 
energy wo, with g being the dimensionless EP coupling 
strength. The Ska (for fc = 1, . . . , N) are the energies of 
noninteracting electrons in the left and right lead (a = 
L,R) with the equilibrium chemical potential /i. The 
corresponding operators c^^ (c^a) create (annihilate) free 
fermions in the N lead states. The last term in Eq. ([1]) 
allows for dot-lead particle transfer. 

We apply to the model (H]) a variational Lang-Firsov 
transformation /^1^^1^^1^^ introducing two parameters 7 
and 7: 

H = sUl)Slh)HSi{^)S2{i) , (2) 
^i(7) = exp{75(6t_fe)dtd}, (3) 
52(7) = exp{7g(&t-6)}. (4) 

S'i(7) describes the antiadiabatic limit, where the 
phononic time scale is much faster than the electronic 
time scale and the deformation of the dot adjusts instan- 
taneously to the presence of an electron. For 7 = 1 it 
coincides with the shifttransformation of the Lang-Firsov 
small-polaron theory, which eliminates the second term 
on the right-hand side of Eq. ([T]) and lowers the dot level 
by the polaron binding energy 

Ep = g^ojQ. (5) 

To account for the competition between polaron local- 
ization and charge transport, an incomplete Lang-Firsov 
transformation with 7 e [0, 1] is used, where 7 will be de- 
termined variationally. The second shift transformation 
S'2(7) describes the regime of fast electron motion, where 
the quasistatic displacement of the equilibrium position 
of the oscillator affects transport. According to similar 
considerations in Ref. [27l the parameter 7 is fixed by 
the condition that the oscillator shift is stationary in the 
equilibrium and steady state. Then 7 = (1 — 7)71^, with 
the dot occupation 

rid = {dU) , (6) 

where (■ ■ • ) denotes the steady state mean value. 
After the transformation the Hamiltonian reads 

H = ?jdU ~ C^{dU - Ud) + Loob^'b + £p(l - jfnl 

k,a k,a 
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with 

A-/i-ep7(2-7)-2ep(l-7)2nd , (8) 

9^19 , ^ka = Ska - fJ. , (9) 

Cka = % e-^(^'-^) , Cd = gcooil - i){b^ + b) . (10) 
V A* 



the external disturbance is exphcit in the time-ordered 
exponential operator 5: 



Here rj is the renormalized energy of the single dot level. 
Cka and Cd are the renormalized interaction coefficients 
of the dot-lead transfer and the EP interaction, respec- 
tively. Note that now the operators d and b represent 
dressed electrons (in analogy to polarons) and the shifted 
local oscillator. The original electron and oscillator op- 
erators, now denoted by d and 6, read 

d = e9(''^-'')d , b = b + gdU+{l--i)gnd. (11) 

We describe the application of a potential difference be- 
tween the leads by adding to Eq. ([7]) the interaction with 
the external fields {/7} and define the voltage bias $ ac- 
cordingly: 

i?i„t =^C/a^cLc,,, with Ua^~5fia. (12) 



$ = (C/l - C/fl)/e , 



(13) 



where e is the (negative) elementary charge. The re- 
sponse of the quantum dot is given by the polaronic 
nonequilibrium real-time Green functions 

9dd{tiM; U) = -i{T du{ti)d\j{t2)) , (14) 
g<^{hM:U) ^i{dl{t2)du{h)) , (15) 
9UhM;U) ^ -i{du{h)dl{t2)) . (16) 

Remember that (■ • • ) denotes the equilibrium average 
with respect to H, while the time dependence of the op- 
erators is now given by _ff -I- Hmt ■ The time ordering 
operator in Eq. (fT4|) is defined by 



Tdu{ti)dlj{t2)^du{ti)dlj{t2) , ti-t2>0 (17) 
= -dlj{t2)du{ti) , ti-i2<0 (18) 

According to Kadanoff-Baym,^^ the real-time response 
functions (fT4l) - (fT6|) may be deduced using the equations 
of motion for the nonequilibrium Green functions of the 
complex time variables t ~ — it, t [0, /3], defined as 

Gdd{ti,t2;U,to) = --^{Trd{h)S{t2)S) , (19) 

GUhM-.UM) - ^^{Trd\t2)d{h)S) , (20) 

G>S,M:UM) = -^{Trd{h)d\t2)S) , (21) 

where the order of ii and t2 is fixed in G^^ and G^^. The 
time dependence of all operators is determined by H and 



(22) 



In Eqs. p9|) -(l2T |) and (|22|) the operator %- orders times 
according to 

Trdu{ti)dl{t2) = d{ti)d\t2) , i(ii - t2) > (23) 
^~S{t2)d{h), i{h-t2)<Q (24) 

In the following, the Green functions of "mixed" opera- 
tors, Gcdik,a;ti,t2;U,to) and gcd{k,a;ti,t2;U), wiU be 
used, which are defined similar to Eqs. (ITil) - (PT|) . The 
functions g follow from the functions G through the lim- 
iting procedure io ^ — oo. 



B. Equations of motion 

We consider the polaronic dot Green function (|19l) . 
where the index "c?d" will be omitted for the moment, 
and start from the Dyson equation in the matrix form 



G(o)-i(ti,t;C/,io)-S(ti,t;i7,io)] •G(f,t2;C/,to) (25) 

^5{tl~t2) . 



In Eq. (|25p the matrix multiplication "•" is defined by 
dt ■ ■ ■ and the S function of complex arguments is 
understood with respect to this integration. With the 
inverse zeroth-order Green function 



G(°)-i(ti,t2) 



i^^~v]Sit,-t2) 



Eq. (HSl) gives for i{ti - to) < i(<2 - ^o) 



i-^^-v)G<{h,t2;U,to) 



(26) 



(27) 



dt I]> {ti,t; U, to)G< {i, t2; U, h) 



to 



+ / dtJ:<{ti,t;U,to)G<{lt2;U,to) 
Jti 

dt E< {ti , t; U, tQ)G> {t, t2;U,to) , 

't2 

where the self-energy functions are defined analo- 
gously to G^: 

S>(ti,i2;t/,to) = S(<i,^2;C/,to), i(ti-i2)>0, (28) 
Y.<{ti,t2;U,to) = T.{ti,t2]U,to), i(ti-i2)<0. (29) 
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On the other hand, the matrix-transposed form of ([25)1 
yields 

[-i^^-v)G<{t,M;UM)^ (30) 
r*i _ 

dt G> {ti,i; U, to)I]< {t, t2; U, to) 

'to 

+ / di G<{tut;U,to)^<it,t2;U,to) 
+ / dtG<lti,t-UMW(tM]UM) ■ 

Jt2 

Similarly to Eqs. (P7|) and (150]) . equations having 
G^{ti,t2; U, to) on the left-hand side are obtained in the 
case i{ti — to) > i(t2 — ^o)- After the limiting procedure 
to — > — oo, we arrive at the equations for the real-time 
response functions of the dot operators: 



i-^^-V^g^ih,t2;U)^ /' [E>(ti,t;C/)-E<(ti;t;t/)].9^(t,t2;;7) (31) 

- r dt J:^{ti,t;U)[g>it,t2;U)-g<{t,t2;U)] , 

^ — OO 

(^~i-^^~V^9^iti,t2;U)= /' [ff>(ti,i;f/)-5<(ii;t;C/)]E^(i,t2;?7) (32) 

dt g^{ti,t;U) [E>{t,t2;U)-i:<{t,t2;U)] . 
The latter equations are general; up to this point no special assumptions or approximations were made. 

C. Steady-state solution 

Limiting ourselves to the steady-state regime, all functions of (^1,^2) will be supposed to depend only ont — ti—t2- 
Then, after suitable change of the integration variables, the difference of the equations for g< in Eqs. ([3T|) and ([32]) 
gives 

CxD 

J dt [g< {i; U)Y.> {t - t; U) - g> {t; U)T.< {t - t;U)] = , (33) 

— 00 

while the differential equation for (g^ —9^) following from Eq. ([?T|) reads 

(i^ - [.9>(i; U) - g<{t; U)] = dt [S>(t; U) - S<(f; U)] [g> {t - t; U) -g<{t- i; U)] (34) 

dt [I]> {t - t; U) ~^<{t- t- U)] [g> {t- U) - g<{t- U)] . 
Using the Fourier transformations of and with factors according to Kadanoff-Baym^^ e.g. 

/oo 
dt g${t;U)e'"' , (35) 
-OO 

/OO 1 
^g^(^;C/)e-'-*, (36) 
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the following exact equations for the steady-state are obtained: 

g<{Lo- U)^>iuj; U) - g> {lo- U)^<{u- U) ^ , 



duj' T,>{uj';U) + Y;<{uj';U) 



[g>{Lu;U)+g<iLu;U)] = 



(37) 
(38) 



OJ — UJ 



If we define, in analogy to the equilibrium expressions 



23 



A{uj;U)^g>{uj;U)+g<{iu;U) , 

f dc. A{u;U) 

9{z;U) = / , 

J it: z — lo 

T{uj-U) = E>(w;C/) + S<(w;f/) , 

dw T{uj]U) 

27r z — UJ ' 



^z-U)-- 
Eq. ([38]) takes the form 



[uj -T^ -ReT.{uj-U)]A{uj;U) 

= T{u- U) Re giuj; U) 

According to Eq. ([M)) we can write 

g<(c^; (7)^^(0.; [/)/>;[/) , 
g>(c^; (7) [/)(!- />;{/)) , 



(39) 
(40) 
(41) 
(42) 

(43) 



(44) 
(45) 



introducing the nonequilibrium distribution /, which fol- 
lows from the steady-state equation ([57]) and the defini- 
tion (gl]) as 



/>; u) 



(46) 



Looking for a solution A{uj\ U) of Eq. (j43|) . which would 
be equal to the equilibrium spectral function for {U} — > 
0, we assume (according to similar considerations in Ref. 
231 ) that g{z;U) has the form 



z — rj — E(z; U) 



(47) 



Together with Eq. gO]), Eq. ^ fulfils Eq. (j43]) identi- 
cally, and the polaronic nonequilibrium spectral function 
becomes 



A{uj;U) 



Tiu;U) 



27T U}—UJ' 



1 2 

+ 


'T(ui;U)' 




2 



(48) 



D. Self-energy 



We determine the polaron self-energy E^d from the 
equations of motion for the generalized Green functions 



of complex time, which were considered for the equilib- 
rium case in Ref. 25. In particular, the coupled equations 
for Gdd and Gcd read 



g'^} \h,t) . G,,(<, h; U, to) = Sih - t2) 



(49) 



+ ^{rrC,{h)d{h)dHt2)S) 



G^^}-\k, a- ti,i; U) . G,^(fc, a; i, U, to) = 

-^{TrCl{tl)d{t,)dHt2)S) , 



(50) 



where, in analogy to Eq. ((26)) . 

Gf^-\k,a-MM;U) = (i^ -4a " u)j S{t,~t2) . 

(51) 

To deduce the functional differential equations for the 

self-energy E^d ~ G^m^^ — G^l ' addition to the phys- 
ical fields {[/}, we introduce the fictitious fields {V} by 
adding to H;^t (cf. Refs. 23, 25, and 28) 

E [^fea(t)C,,(i) + %a{t)CUt)\ + Vd(t)Cd{t) . (52) 

In the same way as in Ref. [2^, the averages on the right- 
hand side of and ([50]) are expressed by means of the 
functional derivatives of Green functions with respect to 
{y}. The resulting functional differential equation for 
Yjdd is solved by iteration to the second order in the in- 
teraction coefficients defined in Eq. pUj) . The correla- 
tion functions of the interaction coefficients are evaluated 
supposing independent Einstein oscillators. Letting then 
{V} — > 0, the following self-consistent result is obtained: 

^dd{tiM:UM) = T.^dliiiM:UM) (53) 
+ [guo{l - i)f Gdd{tiM;UM)Fi{tiM) ■ 
The result of the first iteration step, 

E« (i 1 , ^2 ; {/, to) = E I (^fea > I {k,a-MM;U) (54) 

k.a 

+ Y.\^C^^)\^G'^°Xk,a-MM;U)F,{t,M), 
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is independent of Gdd- The quasiequilibrium nonper- 
turbed Green functions of the leads read 



G(°)<(fc,a;ii,t2;t/) =ie-'«-(*i-*=)/(Cfca + [/a) , (55) 



with f{x) = (e*^^ + 1)-^ The functions Fi and F3 are 
given by and for i(ti — t2) < 0, and by F^ and 
F^ for i(ti — > 0, respectively: 



F' 



fihM) = exp{5^[(nB(c.o) + l)cf-^^o(t,-t,) (gg) 



+ nB(wo)e 



,±iwo(ti-*2) 



(57) 



with nsix) — (e^^ — 1)^^. In Eq. ([5^ . we perform the 
limit io ~^ ~c« and the continuation of the complex time 
variables to real times, while keeping the condition i{ti — 



t2) < for S^^ and i{ti - ^2) > for E^^. We arrive at 



sf,(ti,i2;t/) = S«^(t;{/) 



(58) 



+ [(1- 7)3^0]' .9fd(ti,i2;f/) 
X (ns(wo) + l)e^'"°(*^"*^) 

+ nB(c^o)e^'"°(*^-*^^] , 

^^al^ih M:U) ^Y.\^Cka)? 9i°J^{k, a-MM.U) (59) 

k.a 
X 



where 



{/o(k) +^/s(K)2sinh(s6l) 

S>1 

X (ns(stJo) + l)e=^'''^''^*^-*^) 
+ ni3(swo)e^'^"»(*^-*^)]} , 



^^^'"^^ ^ m!(s + m)! (2) 



(60) 
(61) 



and 



27ri 



A(a;; t/)/(a;; t/) e-'"^*^-*^) 



(62) 



(63) 

Now we insert |(Cfea)P = {\tka\'^ /N) expl-g^ coth0} in 
Eq. (j59p and go from the fc-summation to the integration 
over the lead states with the help of the density of states 
of lead a: 

1 

-El^'^'^l'-'-^E / du^\ta{u^)\'Qa{u^)--- , (64) 

■^^ J. ^ ^ J -00 



(65) 



We then Fourier transform Eq. ((58)) according to Eq. ([35)) 
and, after evaluating the resulting delta functions, obtain 
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(66) 



]«<(a.;L/)=e 



+ [(1 - -/)gujof [a{lu - Wo; U)f{uj - wq; f/)nB(wo) + A{lu + loq; U)f{uj + wq; [/)(ns(wo) + 1) , 

J2 {^o(«)ri°)(^ + + Ua)+J2 sinh(s0) (67) 

a s>l 

{ujo)rf\uj - sujo + ii)f{uj - sujo + Ua) + {neiujo) + l)V'^°\io + sujq + n)J{uj + SUJq + Ua) } , 



-g^ cothe ' 



(68) 



The function I]^^(w; U) can be understood as a generalized in-scattering function of polaron-like quasiparticles at the 
doti^ The second hne in Eq. (p7)) accounts for multiple-phonon emission and, if T > 0, absorption processes. After 
some algebraic manipulations of the Bose- and Fermi- functions, the first-order self-energy (|67p may be written in the 
following form: 



.(i)</ 



(l.; U) = r^l'iio; U)fiio + Ul) + V'^' [u:- U)f{u + Ur) , 
rW(c.;[/) =e-?'™*'''^{ /o(«;)r(")(w + /i)+^/,(.^)2sinh(^ 



(69) 
(70) 



S>1 



:^i°)(w + At - sujo) (riBisuJo) + 1 - f{uj + Ua- swo)) + ri°)(w + /i + sujq) (riBisuJo) + f{uJ + Ua + SCJo)) } 



Because ^^^{ijJ] U) results from interchanging ub ^ (ris + 1), / o (1 - /) and / O (1 - /) in Eqs. (I66l)-(l69]), Eq. (|4T|) 
gives 

r(a;;t/) = rW(c.;[/) (71) 
+ [(1 - i)9UJaY A{uj -ujo;U) (ng (wq) + I - /{uj - ujq; U)^ + A{uj +ujo;U) (ngiujo) + /(w + wq; U)] 



(72) 



From Eq. (|7ip . the spectral function follows using 
Eq. For any parameter 7 < 1, the spectral func- 

tion A and distribution / have to be determined self- 
consistently. Furthermore, because the renormalized dot 
level defined in Eq. ([5]) depends on the dot occupation 
Ud, the latter has to fulfill the self-consistency condition 

/OO 1 
^ />; U)A{lo; U) . (73) 

We note that for 7 = 0, our results are equivalent to the 
SCBA.— For 7 = 1 no self-consistency condition has to 
be fulfilled, as — ^-i^d 77 is independent of n^. 

E. Variational procedure 

To determine the variational parameter 7, we minimize 
the thermodynamic potential which is given by the 
partition function Q as 

17 = -llnQ. (74) 



We assume the leads to be macroscopic objects which 
are negligibly influenced by the states of the dot. Ac- 
cordingly, the contributions of the leads to Q, and to the 
mean energy {H) give only additive constants. Since the 
electronic degrees of freedom of the dot arc coupled to 
the oscillator ones by the second term on the right-hand 
side of Eq. ([7]) , a decoupling approximation will be used 
to determine the electronic part of the thermodynamic 
potential. 

As a consequence of the equation of motion, the fol- 
lowing identity holds: 

?ld\h)d{h) - Cdd\h)d{h) + H'ih) . 

Here H' represents the part of the Hamiltonian ([7]) that 
depends on the operators , d. As an approximation, we 
neglect the second term on the right-hand side of Eq. ([75)) 
and in H'. Taking the statistical averages on both sides 
of Eq. (|75p. remembering that 

{dHh)d{h))^-ig<^{h,t2;U) (76) 
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and using Eq. (|62 



nation of namely 



(H') 



2^ 



(2w - ?j) A{u:; (7)/>; U) (77) 



n' ^n'{X = l) = -^\nQ{X = l) , (79) 
lnQ(A = l) = lnQ(A = 0)-/3^ dXj{V^)^, (80) 



where 



is obtained. To determine the corresponding electronic 
part of the thermodynamic potential, fl' , we consider the 
canonical ensemble given by the Hamiltonian i/^ = Hq + 
V\, where Ho = rjSd and Vx represents the interaction 
part of the Hamiltonian in ([T]) with coefficients XCka and 
XCd, for A e [0, 1]. Applying the result (f77|) gives 



do; 
2^ 



{u:-lT)Ax{u:-U)f{Lo;U) 



(78) 



Here (■ • • )a denotes the dependence of the statistical av- 
erage on A and the indices A on the right-hand side of 
Eq. (I75|) refer to the interaction coefficients in i?^. We 
use the well-known general relations^^''^*' for the determi- 

I 



lnQ(A = 0) = ln(l 



-'7/3 



)■ 



U) 



To make the integration in Eq. (|80|) feasible, the gen- 
eral procedure leading to the thermodynamic poten- 
tial outlined above will be carried out using the so- 
lution for the dot response in the first iteration step, 
described in the preceding section. In particular, the 
spectral function A\{uj;U) is determined according to 
Eq. P5|) . using T^^\uj]U), which is proportional to A^: 
T^^\uj]U) = A2r(i)(w;t7). Similarly, /(w; /7) is deter- 
mined by Eq. using S^j^-*^ and F^^^ on the right-hand 
side. Note however, that will be determined from the 
electron density nd corresponding to the complete self- 
energy I]|;(a;;C/). 

To complete the function J7 which is to be varied with 
respect to 7, we have to take into account the renormal- 
ization of the oscillator energy given in the first line of 
Eq. (O. We finally obtain that 



/3 



/3 



ln(l+e-'"^)+ep(l-7)^n^ + 



ln(l+e-^'5)-|-ep(l-7)^n^ 



dA 



dw 



/(i)(c.)< 



(w-5?)/(i)(w;[/) A2r(i)(w;[/) 



\2'p r duy_ rW{uj';U) 
J 27r bj—uj' 



arctan 



UJ—UJ' 



(82) 



r(i)(tu)/2 



r 



The parameter 7 resulting from the variation of Eq. ([82]) 
is used to determine Y,^j^{uj;U) according to Eq. ([M]). 
The self-energy functions obtained in this way give the 
distribution function f{uj;U) and the spectral function 
A{uj; U) according to Eqs. pS)) and (gS]), respectively. 



of complex times: 



GM{tiM;UM) = -jg^{rrd{h)d\t2)S) (83) 



gib"^ -b){ti) ^~g(b'< ^b)(t2)\ 



F. Relation between electronic and polaronic 
functions 

In the previous sections, the functions A{uj]U) and 
gddi^'j i'^ polaron representation were deduced. Be- 
cause the current through the quantum dot will be given 
by the corresponding electronic functions A{u}; U) and 
'g'^J^uj; U), we have to find a relation between these quan- 
tities. We start by decoupling the fermionic and bosonic 
degrees of freedom in the electronic dot Green function 



Assuming an independent Einstein oscillator, we find 
(7;e^(^'-^)(*i)e-5(^'-'')(*^)) =e-s^'™*''''{/o(At) (84) 



S>1 



where the upper signs correspond to 1(^1 — ^2) > and the 
lower ones to i{ti — ^2) < 0. Going from the complex time 
variables to the real ones, the following relation between 
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U) and g^Jyio; U) is obtained: 

+ 5]/.(K)(e^^gf>±sa;o;f/) 



(85) 



+ e~''5f>Tsc.o;f/))} . 



(86) 
(87) 



With the identities 

e"^ = 2sinh(s6')[l +nB(sa;o)] , 
e"'*^ — 2smh{s6)nB{suJo) , 

the electronic function gf^{uj; U) may be expanded as 

5li(^;C/)=e-^'™*'^^{/o(^)gf,(^;C/) 

+ Is{k)2 sinh(s6l)([l + ns(swo)].9|z(^ ± «^o; C/) 

S>1 



Considering Eqs. p4)) and (|45|) . the electronic spectral 
function is obtained in terms of the polaronic one as 

I(a;; U) = g< (^; {/) + g> (a.; [/) = (89) 
e-s"' '=°*''''|/o(K)A(a;; [/) + ^ /,(k)2 sinh(s6l) 



X [ni3(sa;o) + f{uj + su}q; U)] A{uj + scjq; U) 

+ [nsisujo) + 1 - /(w - swo; ?/)] A{lo - swq; L/)^ | 



G. Current 

The operator of the electron current from lead a to the 
dot reads 



J a = 



(90) 



To calculate the mean value Ja = {Ja), the following con- 
nection of the expectation values to the real-time Green 
functions is used: 



/OO 
— g<i{k,a;uj;U) , 

/OO 1 
-^^[^<(fc,a;a.;[/)] 



(91) 



(92) 



We start from the nonequilibrium Green function of 
the complex time variables for the electron operators, 
namely 

G,d{k,a-hM;U,to)^~^{TrC^^[h)d\t2)S) , (93) 



where S is given by Eq. (|22)) . From the commutators 
with the Hamiltonian in the electron representation, the 
equation of motion is obtained: 

- ika - u}j Gcd{k, a; h,t2;U,to) = (94) 

= Gdd{ii,t2] U, to) . 



(95) 



Equation (|94|) can be rewritten as 

Gcdik,a;ti,t2; U,ta) = 

-% /" " dtG^°J{k,a;h,i;U)Gdd{t,t2;U,to) . 
VN J to 

Performing the limit to —oo while keeping i{ti — 12) < 
0, the following equation for the real-time response func- 
tions is obtained: 



9cdi^.O.\ti,t2;U) 



(96) 



) 

+ I dig(°}<{k,a;h,i;U)g<^it,t2;U) 
dtgi'^:><{k,a;ti,t;U)g<^{t,t2;U) 
di gi°J<{k, a; h,t; U)g>^{t, t2;U) , 



where the quasiequilibrium functions of the noninteract- 
ing leads, gcc'^, coincide with the expressions (|55)) . with 
{ti — ^2) real. Based on Eq. (jHH), the formal manipu- 
lations presented in the Appendix, which are analogous 
to the considerations made in Ref. H, finally lead to the 
following formula for the electron current from the lead 
a to the dot: 

/>oo 

k 

X {fiika + Ua)Ai^; U) - gU^- U)] 

X [f{u + Ua)A{^-,U)~gU^-U)] , 

where the electronic functions g^^{uj;U), A{Ld;U) are 
given by Eqs. (1551) and (15^ . respectively. Since Jl — 
—Jr in steady state, the current formula acquires the 
well-known form^ 



J=1{Jl-Jr) 



(98) 



2 ./-OO 2^ 



r(")(c. + /i) [fLiLu)~fRiLu)]Aiu;U), 
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with fa{i^) = f{uJ + Ua)- In Eq. (l98l) . identical leads 
are assumed, so that r(o)(w) = T^^\uj) = T^^\uj). As a 
check of our numerics, we find indeed that the condition 
Jl — — Jr. holds, as expected for the SCBA. For vanish- 
ing voltage bias $ — > 0, we can express the current as 
J — — L$, where the linear conductance 

L = lim {- J/$} (99) 

<i>-i.O 

results from Eq. ()98|) as 

i - y / 1^ r(°)(.. + m) [-/'(^)] A{u) (100) 

= y'^^rW(c. + ^)/(c.)(l-/(c.))I(c.) 

and the electronic spectral function is now calculated in 
equilibrium. 

III. NUMERICAL RESULTS 

As stated above, the spectral function, dot occupation 
and 7 have to be evaluated self-consistently. We do this 
in a two step manner: (i) for fixed 7 and a starting value 
Ud in Eq. dH) we calculate s|i^j(a;), T'-^\u:). The corre- 
sponding A(^^(cj) and are inserted for A and / in 
the right-hand side of Eqs. and ([7T|) . All functions 
are then iterated until convergence, which is signalled by 

max{\A+iiuj;U) - A,(w; U)\} < 6 , (101) 

with S being a predefined tolerance. In analogy to the 
occurrence of multiple stable solutions in the mean-field 
ansatz of Galperin et al.,^^ for strong EP coupling or high 
voltages, several roots of Eq. ([73)) may exist. We choose 
the root that minimizes the thermodynamic potential, 
(ii) We do this for all parameters 7 to find the global 
minimum of $1(7,71^(7)). The corresponding parameter 
will be referred to as 7min- 

In the following numerical calculations, we suppose 
identical leads and work in the wide-band limit, so that 
is energy independent. 

The equilibrium state, as well as the transport proper- 
ties of molecular junctions crucially depend on the time 
scales of the electronic and phononic subsystem. While 
the lifetime of an electron on the dot is given by the dot- 
lead coupling parameter, t^i oc l/F^^',— the phononic 
time scale is given by the phonon energy Tph oc I/wq- The 
ratio r^^^wo determines which subsystem is the faster 
one. Moreover, one should compare the polaron forma- 
tion time Tpoi oc l/cp to the electron lifetime. If the latter 
is long enough, i.e. if the ratio Sp/T^^^ is large, a tran- 
sient polaron can form at the dot. The parameter will 
yield the mean number of phonons it contains. 



A. Equilibrium situation, low temperature 

We first consider the equilibrium low-temperature 
limit with /x^ = fiji = fi^q = and T = 0.01. Before 
we study the physically more interesting regime of equal 
electronic and phononic time scales, we analyze the two 
limiting cases F^^^ ^ wq and F*^"^ <C ujq. In the following, 
luq = I fixes the energy unit. 

1. Limiting cases 

In the adiabatic case F^°) ^ ujq, the dot deformation 
adjusts quasistatically to the average electronic occupa- 
tion. For small EP coupling, standard perturbation ap- 
proaches are applicable and the expansion of the self- 
energy to second order leads to the Born-approximation 
(BA). On a higher level, the SCBA^^ provides a partial 
resummation of the perturbation series by replacing the 
zero-order Green-function in the BA self-energy with the 
full Green function in a self-consistent way. As was men- 
tioned above, our result ([5^ reduces to the SCBA for 
7^0. 

Figure [IJa) shows the electronic spectral function of 
the adiabatic quantum dot system with A = and Sp — 
5. We compare the SCBA result (7 = 0) to the result 
of the variational calculation, yielding 7i„in = 0.28. The 
SCBA spectrum consists of a single band, whose width 
is given by F^°). Due to the mean-field shift oc = 0.7, 
the renormalized dot level lies beneath the Fermi level 
of the leads (at w = 0) and the dot acts as a tunneling 
well. Because of the short residence time of electrons, 
the effects of inelastic scattering at the dot are small. At 
uj = —loq (lu = +LJo) we find a small peak (dip) in A 
(see arrows) due to narrow logarithmic singularities in 
the denominator of Eq. p8|) .^ 

The variational calculation introduces several correc- 
tions to the spectrum. The finite 7min reduces the ef- 
fective mean-field coupling, i.e. the last term in the po- 
laron shift ([8]). Because it is not fully compensated by 
the rid-independent contribution to Eq. ([8|), the overall 
band shifts upward. In addition, situated at integer mul- 
tiples of Wo from the lead chemical potential, several in- 
elastic resonances form overlapping phononic sidebands. 
Because A{uj = 0) is lowered, transport through the dot 
remains coherent, but with a slightly reduced tunneling 
amplitude. 

In the strong coupling, antiadiabatic case F*^"^ <C 
ujQ, the electron occupies the dot long enough to loose 
coherence and interact with the phonons. Several 
approaches^^^ handle this regime by applying a com- 
plete Fang-Firsov transformation (7 = 1)^ to the Hamil- 
tonian, which gives the exact solution for the isolated 
molecule or when the finite occupation of the leads is 
neglected)^ Consequently, 7min can be considered a mea- 
sure of the small polaron character of the dot state. 

Again we compare the corresponding limit 7 = 1 to 
the result of the variational calculation while setting 
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FIG. 1. For model parameters T = 0.01, fi ^ and $ = 0. 
Panel (a): Electronic spectral functions for r'^-* = 10, Sp = 5, 
A = with 7 = and 7min = 0.28, respectively. Arrows 
mark the phononic features for 7 = 0. Panel (b): Electronic 
spectral functions for F'"' = 0.1, Ep = 1, A = 1 with 7 = 1 
and 7min = 0.81, respectively. Panel (c): Dot occupation and 
variational parameter as functions of the bare dot level A for 
r(°) = 0.1 and e„ = 0.3. 



A = Ep = 1 (see Fig. [Ub)). In the former case, the dot 
level is renormalized by the polaron binding energy and 
represented by the zero-phonon peak at A = A — £p = 0. 
In addition we find pronounced peaks separated by wq, 
signalling the emission of phonons by incident electrons 
and holes. The spectrum documents the formation of a 
long- living polaron state at the dot, with a mean number 



of phonons given by 5^ = 1. 

For the same parameters, the variational calculation 
yields 7i„in = 0.81 < 1 and we find a somewhat broader 
main peak and less spectral weight in the phonon side- 
bands (g^ — 0.66). Consequently, incoherent hopping 
transport through the dot takes place via an intermedi- 
ate polaron state, whose spectral weight and lifetime are 
smaller than predicted by the complete (7 = 1) Lang- 
Firsov calculation. 

Figure [BJc) finally shows the dot occupation and vari- 
ational parameter as functions of the dot level A in the 
antiadiabatic case F*^"^ = 0.1, but for small EP coupling 
Ep = 0.3. In this regime, we find 7i„in ~ 0.7. This 
is in good quantitative agreement with the result of La 
Magna and Deretzis,'''^ who applied a variational Lang- 
Firsov transformation to an effective electron model (cf. 
Fig. 2(b) in Ref. [TtI ). The above calculations show that, 
although the Lang-Firsov approach provides the correct 
physical mechanism, away from the very strong coupling 
limit, adiabatic corrections may not be neglected. 



2. Intermediate dot-lead coupling regime 

We now investigate the regime of comparable electronic 
and phononic time scales by setting F^^^ — 1. Figure [2] 
presents the results of the equilibrium calculation for zero 
to large EP coupling strengths. Shown here are, as func- 
tions of the bare dot level A: the dot occupation (a), 
the variational parameter 7,nin and the renormalized dot 
level rj (b), the linear conductance L (d). For fixed Ep 
and A, Fig. [2jc) gives the thermodynamic potential as 
a function of 7 while Fig. [2](e) and Fig. (f) display the 
electronic spectral functions at A = Ep. 

For Ep = 0, the self-energy ([55]) is exact (black curves 
in Fig. [2]) and the rigid dot acts as a tunneling barrier. As 
A is lowered and the dot charges continuously, the linear 
conductance increases, reaching a maximum at A = 0, 
where the dot level aligns with the lead chemical poten- 
tials and resonant tunneling is possible. The width of 
the conductance resonance is determined by the electron 
lifetime F^"). 

For finite Ep, the variational parameter 7inin ~ 0.5 and 
grows only slightly at A = Ep. As expected, for equal 
electronic and phononic time scales we are far from the 
weak coupling (7 — 0) and strong coupling (7 = 1) limits. 
As a consequence of the EP coupling, the charging tran- 
sition from rid ~ to rid w 1 shifts to higher A because of 
an overall lowering of the effective tunneling barrier. Due 
to the self-consistent mean- field coupling in Eq. ([5]), the 
transition becomes more rapid and even discontinuous 
for Ep > 5 (signalled by the dotted green lines). Here the 
system switches between two stable solutions of Eg. ([75)1 
in analogy to the strong coupling results of Refs. ^M and 
ITtI . Figure [2jc) shows the thermodynamic potential as 
a function of 7 for Ep = 6 with A slightly below and 
above resonance. For 7 < 0.55, the effective mean-field 
coupling in Eq. ([5]) is so strong, that Eq. ([75)) has two 
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FIG. 2. For model parameters P'"^ = 1, T — 0.01, /i = and $ = 0. Panel (a): Dot occupation as a function of the bare dot 
level for several Sp. Panel (b); Variationally determined 7min and renormalized dot level as functions of the bare dot level A. 
Panel (c): thermodynamic potential as a function of 7 for Ep = 6 and A in the vicinity of the discontinuous transition. Here we 
consider the lower (black solid line) or upper (dashed red line) root of the self-consistency equation for Ud- Panel (d): Linear 
conductance as a function of the bare dot level. Panel (e): Electronic spectral function for Ep = 2 at resonance. Panel (f): 
Electronic spectral functions for Sp = 6 and A slightly above (A = G"*") and below (A = 6~) the discontinuous transition. 



roots. For A < Sp, the global minimum of the thermo- 
dynamic potential, situated at 7 = 0.5, corresponds to 
high rid- As A crosses the resonance, the roots change 
roles and the relevant nd jumps. An adiabatic phase 
transition from Ud — to Ud = ^ was also found for a 
single electron at a vibrating quantum dot.^'*''^^ Rapid 
polaron formation and multistability are considered pos- 
sible mechanisms for strongly nonlinear transport prop- 
erties of molecular junctions such as NDC J^"— 

From Fig. ^a.) we see that, in case of a continuous 
transition, Ud = 0.5 whenever A = Sp. As can be eas- 
ily checked from Eq. ([8]), at this point the renormalized 
dot level resonates with the lead chemical potentials, i.e. 
T] = irrespective of 7min- Figure [^Je) shows the corre- 
sponding electronic spectral function for moderate cou- 
pling Ep = A = 2. Few (g^ — 0.5) broad sidebands sig- 
nal phonon emission by either particles (w > 0) or holes 
(w < 0). The spectrum suggests that transmission re- 
mains coherent, but is governed by the slightly increased 
lifetime of the transient polaron state esc l/r(°), with 
r(0) = 0.6. In case of a discontinuous charging, the dot 
level is shifted instantly across the resonance and there is 
no particle-hole symmetric situation, as is demonstrated 
by the spectral functions near the transition for Sp — 6 
(see Fig. [2jf)). Because = 1.5 > 1, spectral weight is 
shifted from the narrow main peak to multiphonon states, 
reducing the tunneling rate in the off-resonant situation 
considerably (Franck-Condon blockade). 



The effects of the EP coupling on the linear response 
of the quantum dot can be seen in Fig. [2l^d). Due to the 
rapid charging and the growing lifetime of the transient 
polaron the symmetrical conductance resonance shifts 
and narrows. This result coincides with the findings of 
Entin-Wohlmann et al. (Ref. [tI) and contradicts the Sp- 
dependent broadening shown in the work of Mitra et al. 
(Ref.^. Note that in case of a continuous transition the 
maximum value of L is independent of the EP coupling 
strength, because the dependence of L on F*^'^-' cancels 
in the low temperature limiti^'^^ In the strong coupling 
limit, the resonance is skipped and the linear response 
signal lowers. In accordance with Refsi^ii, we find no 
side peaks in the linear conductance at low temperatures. 
This is due to "fioating" side bands^ in the electronic 
spectral functions: for all A the phonon signatures are 
offset by wq below and above the lead Fermi level, as can 
be seen from Fig. [IJf). Consequently they are not re- 
solved in the low temperature linear response. This fact 
is missed by single particle approaches;^ 

B. Equilibrium, high temperature 

In the following, we consider the effect of finite tem- 
peratures on the equilibrium properties of the quantum 
dot. We set F^^^ — 0.3 and Ep = 4, thereby entering the 
strong coupling, nonadiabatic regime. Figure |3] shows the 
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FIG. 3. For model parameters P'"' — 0.3, fi — 0, ^ — 0, Sp = i and several temperatures. Panel (a): Dot occupation as a 
function of the bare dot level. Panel (b): Variationally determined 7min and renormalized dot level as functions of the bare 
dot level. Panel (c): thermodynamic potential as a function of 7 for T — 0.3 and A in the vicinity of the resonance. Here we 
consider the lower (black solid line) or upper (dashed red line) root of the self-consistency equation for Ud- Panel (d): Linear 
conductance as a function of the bare dot level. Inset: renormalized dot-lead coupling. Panels (e) and (f): Electronic spectral 
function A and integrated spectral weight S for T = 0.3 and A = 5 and A — 4.5, respectively. 



same quantities as Fig. [21 but compares the low temper- 
ature result (T = 0.01, black curves) to our findings for 
T = 0.3, which, considering phonon energies in the order 
of 100 meVj^i^i corresponds to room temperature. 

Comparing the low temperature result in Fig. ^a) to 
the one for Sp = 4 in Fig. [2{a) we see that the reduction 
of the bare electron tunneling rate increases the effec- 
tive EP coupling strength in such a way that the charg- 
ing transition becomes discontinuous. If we increase the 
temperature the transition becomes continuous again. As 
Fig. (Sljc) shows, for T — 0.3 the optimal 7 is situated in 
a region where only a single root of Eq. ((73l) exists (cf. 
Fig. He)). 

Moreover, at high temperatures the Fermi edges of the 
leads soften. Thermally excited lead electrons see a con- 
siderably reduced injection gap so that the charging tran- 
sition becomes wide spread. We know from Sec. IIII A II 
that in the strong coupling antiadiabatic regime at res- 
onance, when phonon emission by electrons and holes is 
possible, the variational parameter 7,nin comes close to 
unity. At finite temperatures T ^ ujq absorption of free 
phonons by incident electrons opens additional inelastic 
transmission channels. Our ansatz accounts for this with 
7min approaching one at A « 4.5 well above resonance. 
The polaron formation is signalled by two wiggles in the 
renormalized dot level. The impact on the linear con- 
ductance can be seen in Fig. [31Jd): in contrast to the low 
temperature result, we now find three peaks in L. 



Figures |3^e) and (f) compare the electronic spectral 
functions before and after the polaron formation. For 
A = 5 and 7min ~ 0.6, nearly all spectral weight lies 
in a few overlapping emission signals situated above the 
chemical potential. Because at T « wq the floating con- 
dition mentioned in Sec. IIII A^ is relaxed, we find a small 
phonon peak at the chemical potential. That is why the 
conductance resonance broadens with respect to the low 
temperature result. For A — > 4.5, the phonon peaks 
are shifted away from the chemical potential. As 7 ap- 
proaches one, the polaron life time oc l/f^^^ is increased 
by one order of magnitude (see inset Fig. [3Jd)). Conse- 
quently, the peaks in the spectral function narrow and 
spectral weight is transfered to higher order phonon sig- 
nals. The net linear response, being an average over 
transmission channels near the chemical potential, de- 
creases and shapes the outer conductance peaks. At 
A = Ep = 4 the narrow zero phonon peak crosses the 
usual resonance. We note that the maximum value of L 
is smaller than in the low temperature calculation. 



C. Nonequilibrium situation 

The most important experimental technique for the 
characterization of molecular junctions is lETS. Experi- 
ments can be subdivided into nonresonant and resonant 
tunneling scenarios (RIETS). In the former, the energy 
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FIG. 4. For model parameters T = 0.01, F*"' = 10, £p = 2 
and A = 8. Panel (a): Second derivative of the electron 
current as a function of the voltage bias for fixed 7 = (scaled 
by a factor of 20) and variationally determined parameter 
7min, respectively. Panels (b) and (c): Electronic spectral 
functions and their first derivatives at $ — loo- 



of the molecular ion (i.e. rf) lies far above the lead chemi- 
cal potentials. Consequently, electron residence times are 
short and inelastic effects are small. In the latter, reso- 
nance is achieved via the application of a gate voltage 
and strong EP interaction is expected. In both cases 
the current-voltage characteristics exhibit distinct fea- 
tures attributed to vibrational coupling at the junction. 
In analogy to the preceding sections, we will analyze the 
adiabatic and antiadiabatic limiting cases before consid- 
ering equal phononic and electronic time scales. 



1. Limiting cases 

Figure HJ^a) shows the second derivative of the total 
electron current as a function of the voltage in the non- 
resonant (A = 8) adiabatic regime (r(o) = 10) for inter- 
mediate EP coupling strength {sp — 2). For fixed 7 = 
we find a single dip at $ = wqj where rj = 6.8. Here, 
phonon emission by incident electrons causes an addi- 
tional inelastic tunneling current. Moreover, quasielastic 
processes involving the emission and subsequent absorp- 



tion of a single phonon are no longer virtual, because 
the intermediate polaron state is only partially occupied. 
The tunneling current (|98p is an integral over the ener- 
gies of all incident and outgoing electrons and does not 
resolve the various tunneling processes. Therefore pola- 
ronic features are observed in the second derivative of 
J. As Persson showed;^ the destructive interference of 
the elastic and quasielastic processes may overcompen- 
sate the positive inelastic contribution, leading to the 
dip in the lETS signal. In their SCBA analysis, Galperin 
et al. (Ref. 11) demonstrated the strong qualitative de- 
pendence of this signature on the dot level A and the 
bare molecule-lead coupling T^^\ Our ansatz allows for 
the polaronic renormalization of both these parameters: 
At $ = ojo the variational calculation gives an optimal 
7min = 0.3 and the effective dot level is further lowered 
[rj — 6.4 at $ = wq). As can be seen from the electronic 
spectral function in Fig. Hl^b), the spectral weight of in- 
elastic electron tunneling processes at cj > $/2 — 0.5 
grows at the cost of the elastic transmission at w = 0. As 
a consequence, the overall lETS signal now shows a pro- 
nounced peak at $ = wq (note the scaling of the curves in 
Fig. m^a)) and additional phonon features whenever the 
voltage crosses integer multiples of loq. With the current 
being an integral over the quantum dot spectrum, the 
qualitative change in the one-phonon lETS signal can be 
traced back the first derivative of A{lo),^^ which can be 
seen in Fig. SJ^c). When going from 7 = 1 to 7inin — 0.3, 
the sum of the peak derivatives oi AaXuj — ^jll.r — ±<&/2 
changes sign, showing that the inelastic tunneling current 
outweighs the destructive interference of the elastic chan- 
nels. 

Figure a) and (b) present the total current and dif- 
ferential conductance as functions of the voltage in the 
resonant (A = 2) antiadiabatic regime (r(o) = 0.1) for 
intermediate EP coupling strength [sp = 2). Because 
the voltage is raised symmetrically around the equilib- 
rium chemical potential the dot occupation as well as the 
renormalized dot level 77 = remain constant. Both, the 
variational calculation and the 7=1 case exhibit steps in 
the total current and pronounced peaks in the differential 
conductance whenever the voltage equals multiple inte- 
gers of 2u)Q. Here resonant tunneling through phononic 
sidebands becomes possible. At $ « 12, the current 
saturates because now the so-called "Fermi window" 
a; g [—<i>/2, -!-<!>/ 2] encompasses all phonon side bands 
(see. Fig. 5(c)). In the low-voltage region $ < 4, the 
optimal variational parameter differs considerably from 
one (7niin ~ 0.9), thereby increasing the overall weight of 
the relevant few-phonon inelastic tunneling channels. As 
a consequence, the low-voltage current is larger than in 
the 7 = 1 case. Nevertheless, the growth of 7ijiin along a 
current plateau dynamically shifts spectral weight from 
the corresponding resonant inelastic channel to higher ly- 
ing bands outside the Fermi window. As can be seen from 
the inset of Fig. [Hb) , the differential conductance is neg- 
ative, which is in accordance with the polaron induced 
NDC found by La Magna and Deretzis.^^ Only when an 
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FIG. 5. For model parameters T = 0.01, r'°> = 0.1, £p = 2 
and A = 2. Panel (a): Electron current as a function of the 
voltage bias, compared to the result with fixed 7 = 1. Inset: 
renormalized dot-lead coupling. Panel (b): Differential con- 
ductance as a function of the voltage bias. Inset: Zoom on 
the low- voltage region. Panel (c): Electronic spectral func- 
tions A and nonequilibrium electron distribution functions / 
for several voltages. 
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set: Zoom on the low- voltage region. Panel (c): Electronic 
spectral functions A and nonequilibrium electron distribution 
functions / for several voltages. 
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upward step (peak in d^J/d$^) signals the opening of 
a nonresonant inelastic channel, the differential conduc- 
tance becomes positive again. 



2. Intermediate dot-lead coupling regime 

We now turn to the regime of equal electronic and 
phononic time scales, setting T^^'' — 1 and keeping 
T = 0.01 and Ep = 2 fixed. First, we hold A = 2 at 
resonance, starting with 7,nin = 0.5 and — 0.5 in 
equilibrium (cf. Fig. [2]). Figure [BJa) presents the cor- 
responding current-voltage characteristics. We compare 
the result of the variational calculation (black solid lines) 
to the case with fixed 7 = 1 (blue dashed lines) and to an 
effective electron model (red dash-dotted lines) . The lat- 
ter is obtained by setting g = in Eqs. and ((7T]) and 

inserting for ri"' the renormalized dot- lead coupling r*^"^ 
resulting from the variational calculation. It is compara- 
ble to earlier works where the averaging over the phonon 
state leads to an effective electron HamiltonianJii^ 

With growing voltage the variational parameter 
steadily increases and approaches one in the high- voltage 
limit <J> > 6. The elastic transmission rate r*^"^ shown 
in the inset of Fig. IHl^a) decreases accordingly. It ex- 
hibits steps at integer multiples of 2woj suggesting that 
the polaron formation is especially rapid whenever a new 
resonant inelastic channel is accessible. The electronic 
spectral functions in Fig. [Sfc) show that spectral weight 
is shifted from the zero-phonon peak to the overlapping 
phonon side bands. 

The current-voltage characteristics of the interacting 
results (7min and 7 = 1) contain signatures of both lim- 
iting cases discussed in Sec. IIII C 11 as can be seen from 
the differential conductance in Fig. IH^b). As before, at 
voltages corresponding to integer multiples of 2a;o, steps 
in the current (peaks in the conductance) signal the onset 
of resonant inelastic tunneling. These steps are consid- 
erably broadened and overlap with the onset of nonreso- 
nant inelastic tunneling. As a consequence, the polaron 
induced renormalization of the resonant channel is com- 
pensated and, in contrast to the low- voltage antiadiabatic 
regime, dJ/d^ remains strictly positive. 

The effective electron model overestimates the current 
in the region < $ < 2a;o. Since the spectrum con- 
tains no phonon side bands, for $ > 2a;o the decrease 
of the elastic tunneling rate oc V^^^ is not compensated 
by resonant or nonresonant inelastic transmission pro- 
cesses. Consequently, we find a considerably lower max- 
imum current and, in accordance with the results of La 
Magna and Deretzis^^ NDC in the intermediate-to-high 
voltage region. We conclude that the polaron induced 
renormalization of the dot-lead coupling is indeed a pos- 
sible mechanism for NDC. Yet, the effective electron cal- 
culation misses the spectral features that are essential for 
electron transport at voltages exceeding loq. The inter- 
play of several inelastic transmission channels may heav- 
ily reduce or, for F^^^ > wq, even prevent the occurrence 
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of NDC. 

Another interesting consequence of the dynamic po- 
laron formation can be observed in the high voltage 
regime, where a crossover from nonresonant to resonant 
transport takes place. We keep the above system param- 
eters, but start from the nonresonant equilibrium situa- 
tion with A = 8. The result is presented in Fig. [T] As 
the voltage is raised, the variational parameter as well as 
the effective dot level remain nearly constant and trans- 
port takes place via nonresonant inelastic tunneling. At 
$ — 12.4 the chemical potential of one lead resonates 
with 77 = 6.2, causing a broad step in the total current. 
When the voltage is raised further, the system maximizes 
its kinetic energy by decreasing the polaronic shift in such 
a way, that 77 stays locked to the lead chemical potential 
(see Fig. Efb)). As the spectral functions in Fig. [71[c) 
suggest, this happens at the cost of the inelastic trans- 
mission channels. As soon as 7inin — and resonance of 
the zero-phonon level can no longer be maintained, the 
system reduces its potential energy by forming a tran- 
sient polaron. Here 7min jumps to 1 and the effective dot 
level is lowered by the full polaron binding energy e^. 
The spectral functions in the vicinity of this transition 
show that the spectral weight is redistributed to inelas- 
tic channels within the Fermi window. Consequently, the 
current shows no discontinuity or NDC at this point. 



IV. SUMMARY 

In this work, we investigate the steady-state transport 
through a vibrating molecular quantum dot. Within the 
Kadanoff-Baym formalism, the nonequilibrium dot self- 
energy is calculated to second order in the interaction 
coefficients. To describe the polaronic character of the 
quantum dot state, we apply a variational Lang-Firsov 
transformation and determine the degree of transforma- 
tion self-consistently by minimizing the thermodynamic 
potential. 

In this framework we are able to study the molecular 
junction for all ratios of the dot-lead coupling to the en- 
ergy of the local phonon mode, i.e. from the adiabatic 
to the antiadiabatic regime. Moreover, the EP interac- 
tion can be varied from weak to strong coupling. Tuning 
the electronic dot level and the external voltage bias, we 
can finally consider resonant and off-resonant transport 
in the equilibrium and nonequilibrium situation. 

In the adiabatic regime, we find important correc- 
tions to the result of the SCBA when the EP coupling 
grows: In the equilibrium, off-resonant situation, the 
mean-field oscillator shift is reduced and spectral weight 
is transferred from elastic to inelastic channels. For finite 
voltages, we observe a pronounced peak in the electron 
tunneling signal, followed by several pronounced multi- 
phonon features. 

In the antiadiabatic regime, away from the very strong 
coupling limit, the weight of the transient polaron state 
is smaller than predicted by the complete Lang-Firsov 



transformation. Accordingly, the equilibrium linear con- 
ductance as well as the low voltage resonant tunneling 
current increase, because few-phonon emission processes 
are amplified. As the voltage bias grows the full Lang- 
Firsov polaron forms. Here, due to a dynamical renor- 
malization of the dot-lead coupling, we find NDC along 
the resonant current plateaus. 

Most notably, our variational approach also allows the 
investigation of the intermediate regime where the dot- 
lead coupling and the phonon energy are of the same or- 
der. For weak EP coupling, the linear conductance shows 
a single resonance peak as a function of the electronic dot 
level. When the coupling strength is increased this peak 
narrows and shifts, signaling the crossover from coherent 
tunneling to sequential hopping via a long-living, tran- 
sient polaron at the dot. For very strong coupling, the 
polaron formation takes place discontinuously, as the sys- 
tem switches between various metastable states. At finite 
temperatures, this transition becomes continuous again. 
At the same time, the equilibrium linear conductance 
signal broadens and shows distinct phonon side peaks. 
Thermally activated transport via phonon absorption in- 
duces polaron formation far from resonance. In the low- 
temperature, nonequilibrium situation, the differential 
conductance remains positive for all voltages: the po- 
laron induced renormalization of the dot-lead coupling is 
compensated by the onset of off-resonant inelastic trans- 
port. In the off-resonant, high-voltage regime, the po- 
laron level follows the lead chemical potential to enhance 
resonant transport and maximize the kinetic energy. 

Let us emphasize that we determine the current 
through the dot by means of an approximation to the 
electronic spectral function that contains inelastic fea- 
tures to all orders in the EP coupling. We compare our 
results to an effective electron model, which accounts for 
the electron-phonon interaction only via a renormalized 
dot-lead coupling parameter (e.g. in analogy to Ref. [TtI)- 
For this model negative differential conductance is ob- 
served. This is because the effective electronic spectral 
function does not include inelastic features that affect 
transport for voltages exceeding the phonon frequency. 

The present study may be extended in several di- 
rections: (i) description of hysteretic behavior in the 
strong coupling, high voltage regime; (ii) inclusion of 
the dynamics of the phonon subsystem by means of 
nonequilibrium phonon Green functions; (iii) incorpora- 
tion of Coulomb interaction at the dot to produce even 
stronger nonlinear effects through the competition of a 
population-dependent repulsive dot potential with the 
polaronic level shift. 



ACKNOWLEDGMENTS 

This work was supported by Deutsche Forschungsge- 
meinschaft through SFB 652 B5. TK and HE acknowl- 
edge the hospitality at the Institute of Physics ASCR. 



18 



APPENDIX: DERIVATION OF THE CURRENT 
FORMULA 



Deducing the current response in Sec. Ill G| the fol- 
lowing real-time Green functions (defined according to 
Mahan^i) are used: 



AhM) = Q{t2-ti)g>{ti,t2) 
+ Q{h-t2)g<{ti,t2) 



(102) 
(103) 



where Q is the Heaviside function. The relations of 5* and 
g* to the retarded and advanced Green functions read 



and Eq. (I96p may be written as 



-—g<^{k,a;ti,t2;U) 



(104) 
(105) 



(106) 



di 1 5^'* (fc, a; i 1 , fi ; U)g<^ih M^U) 
dh gi°J<{k,a;h,h;U)gld{ti,t2;U) . 



As far as the steady-state is concerned, all averages in 
the definitions of the Green functions above dependent 

I 



only on the differences of time variables. Consequently, 
the integrals on the right-hand side of Eq. (|106p may be 
rewritten in the form of a convolution and the Fourier 
transformation of Eq. (|106p is 



t 



ka 

N 



£^'{k,a;io;U)g<,iio;U) (107) 



/,'!}<ik,a;uj;U)gUu;U) 



Here, the Fourier transforms of the response functions are 
defined in the usual convention, i.e. without the factors 
±i introduced by Eqs. (|35l) and (|36l) . In particular, the 
conventional Fourier transforms fulfil 



-.9^(^) , 



(108) 



because the left-hand side of Eq. ([35|) is a real func- 
tion. Taking into account the general property that 



[g-*(^)] 
give 



9 



adv 



(w), the relations (fT04| . (fT05|) and (fTOS]) 



[9\^)Y 



(109) 



With the help of Eqs. pII5)) and I^TU^ . the complex con- 
jugate of 5^(fc, a; uj; U) in Eq. is determined and the 
following formula for the current Jq results: 



-g(?<(fc,a;c.;C/)[g*,,(u;;C/) + g*,(^;C/)] } . 



(110) 



Substituting the explicit forms of the free electron functions ^c?^ and using the relation g^ -\- g^ = g^ + g^ following 
from Eqs. (I104p and (jlOSI) . we obtain 



= E I ^ 27rJ(^ - 6a){ - ig^A^; U) + fiCka + Ua)i[g^,{uj; U) - g>,{uj; U)] } . (Ill) 



Going back to the definitions of the Fourier transforms according to Eqs. (1551) and we arrive at Eq. (1^71) of 
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